home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / result.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  8.8 KB  |  281 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module result)
  13.  
  14. (DECLARE-TOP(SPECIAL VARLIST GENVAR $RATFAC $KEEPFLOAT MODULUS *ALPHA XV))
  15.  
  16. (LOAD-MACSYMA-MACROS RATMAC)
  17.  
  18. (DECLARE-TOP(SPLITFILE MRESUL))
  19.  
  20. (DEFMFUN $POLY_DISCRIMINANT (POLY VAR)
  21.   (LET* ((VARLIST (LIST VAR))
  22.      (GENVAR ()) 
  23.      (RFORM (RFORM POLY))
  24.      (RVAR (CAR (LAST GENVAR)))
  25.      (N (PDEGREE (SETQ POLY (CAR RFORM)) RVAR)))
  26.     
  27.     (COND ((= N 1) 1)
  28.       ((OR (= N 0) (NOT (atom  (CDR RFORM))))
  29.        (MERROR "ARG. MUST BE A POLYNOMIAL IN VAR"))
  30.       (T (PDIS (PRESIGN
  31.             (// (f* N (f1- N)) 2)
  32.             (PQUOTIENT (RESULTANT POLY (PDERIVATIVE POLY RVAR))
  33.                    (P-LC POLY))))))))
  34.  
  35. (DEFMFUN $RESULTANT (A B MAINVAR)
  36.        (PROG (VARLIST FORMFLAG $RATFAC RES ANS GENVAR $KEEPFLOAT)
  37.          (SETQ VARLIST (LIST MAINVAR) $RATFAC T ANS 1)
  38.          (AND ($RATP A)(setq formflag t)(SETQ a ($ratdisrep a)))
  39.          (AND ($RATP B)(setq formflag t)(SETQ b ($ratdisrep b)))
  40.          (NEWVAR A)
  41.          (NEWVAR B)
  42.          (SETQ A (LMAKE2 (CADR (RATREP* A)) NIL))
  43.          (SETQ B (LMAKE2 (CADR (RATREP* B)) NIL))
  44.          (SETQ MAINVAR (CAADR (RATREP* MAINVAR)))
  45.        (DO ((L1 A (CDR L1))) ((NULL L1))
  46.        (DO ((L2 B (CDR L2))) ((NULL L2))
  47.            (SETQ RES (RESULT1 (CAAR L1) (CAAR L2) MAINVAR))
  48.            (SETQ ANS (PTIMES ANS (PEXPT 
  49.          (COND ((ZEROP (CADDR RES)) (CAR RES))
  50.                (T (PTIMESCHK (CAR RES)
  51.                      (PEXPT (MAKPROD (CADR RES) NIL)
  52.                         (CADDR RES)))))
  53.          (TIMES (CDAR L1) (CDAR L2)))))))
  54.     (RETURN (COND (FORMFLAG (PDIS* ANS)) (T (PDIS ANS))))))
  55.  
  56. (DEFUN RESULT1 (P1 P2 VAR)
  57.        (COND ((OR (PCOEFP P1) (POINTERGP VAR (CAR P1)))
  58.           (LIST 1 P1 (pdegree P2 VAR)))
  59.          ((OR (PCOEFP P2) (POINTERGP VAR (CAR P2)))
  60.           (LIST 1 P2 (pdegree P1 VAR)))
  61.          ((NULL (CDDDR P1))
  62.           (COND ((NULL (CDDDR P2)) (LIST 0 0 1))
  63.             (T (LIST (PEXPT (CADDR P1) (CADR P2))
  64.                  (PCSUBSTY 0 VAR P2)
  65.                  (CADR P1)))))
  66.          ((NULL (CDDDR P2))
  67.           (LIST (PEXPT (CADDR P2) (CADR P1))
  68.             (PCSUBSTY 0 VAR P1)
  69.             (CADR P2)))
  70.          ((> (SETQ VAR (GCD (PGCDEXPON P1) (PGCDEXPON P2))) 1)
  71.           (LIST 1 (RESULTANT (PEXPON*// P1 VAR NIL)
  72.                  (PEXPON*// P2 VAR NIL)) VAR))
  73.          (T (LIST 1 (RESULTANT P1 P2) 1))))
  74.  
  75. (DEFMVAR $RESULTANT '$SUBRES "Designates which resultant algorithm")
  76.  
  77. (DEFVAR *resultlist '($subres $mod $red))
  78.  
  79. (DEFMFUN RESULTANT (P1 P2)                ;assumes same main var
  80.   (IF (> (P-LE P2) (P-LE P1))
  81.       (PRESIGN (f* (P-LE P1) (P-LE P2)) (RESULTANT P2 P1)) 
  82.       (CASE $RESULTANT
  83.     ($SUBRES (SUBRESULT P1 P2))
  84.     #+broken ($MOD (MODRESULT P1 P2))
  85.     ($RED (REDRESULT P1 P2))
  86.     (T (MERROR "No such resultant algorithm")))))
  87.  
  88. (DEFUN PRESIGN (N P)
  89.   (IF (ODDP N) (PMINUS P) P))
  90.  
  91. (declare-top (SPLITFILE SUBRES))
  92. ;computes resultant using subresultant p.r.s. TOMS Sept. 1978
  93.  
  94. (defun subresult (p q)
  95.   (sloop for g = 1 then (p-lc p)
  96.     for h = 1 then (pquotient (pexpt g d) h^1-d)
  97.     for degq = (pdegree q (p-var p))
  98.     for d = (f- (p-le p) degq)
  99.     for h^1-d = (if (equal h 1) 1 (pexpt h (f1- d)))
  100.     if (zerop degq) return (if (pzerop q) q (pquotient (pexpt q d) h^1-d))
  101.     do (psetq p q
  102.           q (presign (f1+ d) (pquotient (prem p q)
  103.                            (ptimes g (ptimes h h^1-d)))))))
  104.  
  105. (DECLARE-TOP (SPLITFILE REDRES))
  106.  
  107. ;    PACKAGE FOR CALCULATING MULTIVARIATE POLYNOMIAL RESULTANTS 
  108. ;    USING MODIFIED REDUCED P.R.S.
  109.  
  110. (DEFUN REDRESULT (U V)
  111.    (PROG (A R SIGMA C)
  112.      (SETQ A 1)
  113.      (SETQ SIGMA 0)
  114.      (SETQ C 1)
  115.     A    (IF (PZEROP (SETQ R (PREM U V))) (RETURN (PZERO)))
  116.      (SETQ C (PTIMESCHK C (PEXPT (P-LC V)
  117.                      (f* (f- (P-LE U) (P-LE V))
  118.                     (f- (P-LE V) (PDEGREE R (P-VAR U))
  119.                        1)))))
  120.      (SETQ SIGMA (f+ SIGMA (f* (P-LE U) (P-LE V))))
  121.      (IF (ZEROP (PDEGREE R (P-VAR U)))
  122.          (RETURN
  123.           (PRESIGN SIGMA
  124.                (PQUOTIENT (PEXPT (PQUOTIENTCHK R A) (P-LE V)) C))))
  125.      (PSETQ U V
  126.         V (PQUOTIENTCHK R A)
  127.         A (PEXPT (P-LC V) (f+ (P-LE U) 1 (f- (P-LE V)))))
  128.      (GO A)))
  129.  
  130.  
  131. (declare-top (SPLITFILE MODRES))
  132.  
  133. ;    PACKAGE FOR CALCULATING MULTIVARIATE POLYNOMIAL RESULTANTS 
  134. ;    USING MODULAR AND EVALUATION HOMOMORPHISMS.
  135. ; modresultant fails on the following example
  136. ;RESULTANT(((-4)*Z)^4+(Y+8*Z)^4+(X-5*Z)^4-1,
  137. ;           ((-4)*Z)^4-(X-5*Z)^3*((-4)*Z)^3+(Y+8*Z)^3*((-4)*Z)^2
  138. ;             +(-2)*(Y+8*Z)^4+((-4)*Z)^4+1,Z)
  139.  
  140. #+broken
  141. (progn
  142. (DEFUN MODRESULT (A B)
  143.        (MODRESULT1 A B (SORT (UNION* (LISTOVARS A) (LISTOVARS B))
  144.                  (FUNCTION POINTERGP))))
  145.  
  146. (DEFUN MODRESULT1 (X Y VARL)
  147.        (COND ((NULL MODULUS) (PRES X Y (CAR VARL) (CDR VARL)))
  148.          (T (CPRES X Y (CAR VARL) (CDR VARL))) ))
  149.  
  150. (DEFUN PRES (A B XR1 VARL)
  151.        (PROG (M N F A* B* C* P Q C MODULUS HMODULUS)
  152.     (SETQ M (CADR A))
  153.     (SETQ N (CADR B))
  154.     (SETQ F (COEFBOUND M N (MAXNORM (CDR A)) (MAXNORM (CDR B)) ))
  155.     (SETQ Q 1)
  156.     (SETQ C 0)
  157.     (SETQ P *ALPHA)
  158.     (GO STEP3)
  159. STEP2    (SETQ P (NEWPRIME P))
  160. STEP3    (SETQMODULUS P)
  161.     (SETQ A* (PMOD A))
  162.     (SETQ B* (PMOD B))
  163.     (COND ((OR (REJECT A* M XR1) (REJECT B* N XR1)) (GO STEP2)))
  164.     (SETQ C* (CPRES A* B* XR1 VARL))
  165.     (SETQMODULUS NIL)
  166.     (SETQ C (LAGRANGE3 C C* P Q))
  167.     (SETQ Q (TIMES P Q))
  168.     (COND ((GREATERP Q F) (RETURN C))
  169.           (T (GO STEP2)) ) ))
  170.  
  171. (DEFUN REJECT (A M XV)
  172.        (NOT (EQN (PDEGREE A XV) M)))
  173.  
  174. (DEFUN COEFBOUND (M N D E)
  175.        (TIMES 2 (EXPT (f1+ M) (// N 2))
  176.           (EXPT (f1+ N) (// M 2))
  177.           (COND ((ODDP N) (f1+ ($ISQRT (f1+ M))))
  178.             (T 1))
  179.           (COND ((ODDP M) (f1+ ($ISQRT (f1+ N))))
  180.             (T 1))
  181. ; (FACTORIAL (PLUS M N)) USED TO REPLACE PREV. 4 LINES. KNU II P. 375
  182.           (EXPT D N)
  183.           (EXPT E M) ))
  184.  
  185. (DEFUN MAIN2 (A VAR EXP TOT)
  186.        (COND ((NULL A) (CONS EXP TOT))
  187.          (T (MAIN2 (CDDR A) VAR
  188.                (MAX (SETQ VAR (PDEGREE (CADR A) VAR)) EXP)
  189.                (MAX (f+ (CAR A) VAR) TOT))) ))
  190.  
  191. (DEFUN CPRES (A B XR1 VARL)                ;XR1 IS MAIN VAR WHICH
  192.        (COND ((NULL VARL) (CPRES1 (CDR A) (CDR B)))    ;RESULTANT ELIMINATES
  193.          (T    (PROG (  M2           ( M1 (CADR A))
  194.           ( N1 (CADR B))  N2 (K 0) C D A* B* C* BP XV);XV IS INTERPOLATED VAR
  195.               (DECLARE (FIXNUM M1 N1 K))
  196.  
  197.             STEP2
  198.           (SETQ XV (CAR VARL))
  199.           (SETQ VARL (CDR VARL))
  200.           (SETQ M2 (MAIN2 (CDR A) XV 0 0))    ;<XV DEG . TOTAL DEG>
  201.           (SETQ N2 (MAIN2 (CDR B) XV 0 0))
  202.           (COND ((ZEROP (f+ (CAR M2) (CAR N2)))
  203.              (COND ((NULL VARL) (RETURN (CPRES1 (CDR A) (CDR B))))
  204.                    (T (GO STEP2)) ) ))
  205.           (SETQ K (f1+ (MIN (f+ (f* M1 (CAR N2)) (f* N1 (CAR M2)))
  206.                    (f+ (f* M1 (CDR N2)) (f* N1 (CDR M2))
  207.                       (f- (f* M1 N1))) )))
  208.           (SETQ C 0)
  209.           (SETQ D 1)
  210.           (SETQ M2 (CAR M2) N2 (CAR N2))
  211.           (SETQ BP (MINUS 1))
  212.            STEP3
  213.               (COND ((EQUAL (SETQ BP (ADD1 BP)) MODULUS)
  214.              (merror "Resultant primes too small."))
  215.             ((ZEROP M2) (SETQ A* A))
  216.             (T (SETQ A* (PCSUBST A BP XV))
  217.                (COND ((REJECT A* M1 XR1)(GO STEP3)) )) )
  218.           (COND ((ZEROP N2) (SETQ B* B))
  219.             (T (SETQ B* (PCSUBST B BP XV))
  220.                (COND ((REJECT B* N1 XR1) (GO STEP3))) ))
  221.           (SETQ C* (CPRES A* B* XR1 VARL))
  222.           (SETQ C (LAGRANGE33 C C* D BP))
  223.           (SETQ D (PTIMESCHK D (LIST XV 1 1 0 (CMINUS BP))))
  224.           (COND ((> (CADR D) K) (RETURN C))
  225.             (T (GO STEP3))) )) ))
  226. )
  227.  
  228. (declare-top (SPLITFILE BEZOUT))
  229.  
  230. ;; *** NOTE THAT MATRIX PRODUCED IS ALWAYS SYMETRIC
  231. ;; *** ABOUT THE MINOR DIAGONAL.
  232.  
  233. (DEFMFUN $BEZOUT (P Q VAR)
  234.        (LET ((VARLIST (LIST VAR)) GENVAR)
  235.         (NEWVAR P)
  236.         (NEWVAR Q)
  237.         (SETQ P (CADR (RATREP* P))
  238.           Q (CADR (RATREP* Q)))
  239.         (SETQ P (COND ((> (CADR Q) (CADR P)) (BEZOUT Q P))
  240.               (T (BEZOUT P Q))))
  241.         (CONS '($MATRIX)
  242.           (MAPCAR #'(LAMBDA (L) (CONS '(MLIST) (MAPCAR 'PDIS L)))
  243.               P))))
  244.  
  245. (DEFUN VMAKE (POLY N *L)
  246.        (DO ((I (f1- N) (f1- I))) ((MINUSP I))
  247.        (COND ((OR (NULL POLY) (< (CAR POLY) I))
  248.           (SETQ *L (CONS 0 *L)))
  249.          (T (SETQ *L (CONS (CADR POLY) *L))
  250.             (SETQ POLY (CDDR POLY)))))
  251.        (NREVERSE *L))
  252.  
  253. (DEFUN BEZOUT (P Q)
  254.   (LET* ((N (f1+ (P-LE P)))
  255.      (N2 (f- N (P-LE Q)))
  256.      (A (VMAKE (P-TERMS P) N NIL))
  257.      (B (VMAKE (P-TERMS Q) N NIL))
  258.      (AR (REVERSE (NTHCDR N2 A)))
  259.      (BR (REVERSE (NTHCDR N2 B)))
  260.      (L (NZEROS N NIL)))
  261.     (RPLACD (NTHCDR (f1- (P-LE P)) A) NIL)
  262.     (RPLACD (NTHCDR (f1- (P-LE P)) B) NIL)
  263.     (NCONC
  264.      (MAPCAR
  265.       #'(LAMBDA (AR BR)
  266.       (SETQ L (MAPCAR #'(LAMBDA (A B L)
  267.                   (PPLUSCHK L (PDIFFERENCE
  268.                        (PTIMES BR A) (PTIMES AR B))))
  269.               A B (CONS 0 L))))
  270.       AR BR)
  271.      (AND (PZEROP (CAR B))
  272.       (DO ((B (VMAKE (CDR Q) (CADR P) NIL) (ROT* B))
  273.            (M NIL (CONS B M)))
  274.           ((NOT (PZEROP (CAR B))) (CONS B M))))) ))
  275.  
  276. (DEFUN ROT* (B)
  277.   (SETQ B (COPY1 B)) (PROG2 (NCONC B B) (CDR B) (RPLACD B NIL)))
  278.  
  279.  
  280. (DEFUN PPLUSCHK (P Q) (COND ((PZEROP P) Q) (T (PPLUS P Q))))
  281.